%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% File:               bootstrap_game.m
%
% Author:             Miguel R. Rueda and Sergio Ascencio
%
% Description:        Computes boostrap standard errors using 2 step ML estimator. 
% Language:           MATLAB R2013b (8.2.0.701) 64 Bit

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [betas,se,pvals,UCI,LCI]=bootstrap_game(Y,X,beta_p,nReps,state,CLv_PAN,CLv_PRI,CLv)

N=size(X,1);
cols=length(beta_p(:));
betas=zeros(cols,nReps);

if nargin<8
    rng(1e3)
    id = ceil(rand(N,nReps)*N);    
    for i=1:nReps    
        [beta,~,~,~,~,~]=representatives_game(Y(id(:,i),:),X(id(:,i),:),state,CLv_PAN(id(:,i),:),CLv_PRI(id(:,i),:),beta_p);
        betas(:,i)=beta(:);
    end
    
else 
    CLv_u=unique(CLv);
    N_cl=length(CLv_u);
    rng(1e3)
    id_aux = ceil(rand(N_cl,nReps)*N_cl);
    for i=1:nReps
        rn=CLv_u(id_aux(:,i),1);
        rn_u=unique(rn);
        out = [rn_u,histc(rn(:),rn_u)];
        Yb_c=cell(length(out(:,1)),1);
        Xb_c=cell(length(out(:,1)),1);
        CLv_PAN_c=cell(length(out(:,1)),1);
        CLv_PRI_c=cell(length(out(:,1)),1);
        tic;
        for j=1:length(out(:,1))
            Yb_c{j,1}=kron(Y(CLv==out(j,1),:),ones(out(j,2),1));
            Xb_c{j,1}=kron(X(CLv==out(j,1),:),ones(out(j,2),1));
            CLv_PAN_c{j,1}=kron(CLv_PAN(CLv==out(j,1),:),ones(out(j,2),1));
            CLv_PRI_c{j,1}=kron(CLv_PRI(CLv==out(j,1),:),ones(out(j,2),1));
        end 
        Yb=cell2mat(Yb_c);
        Xb=cell2mat(Xb_c);
        CLv_PANb=cell2mat(CLv_PAN_c);
        CLv_PRIb=cell2mat(CLv_PRI_c);
        toc;
        [beta_it,~,~,~,~,~]=representatives_game(Yb,Xb,state,CLv_PANb,CLv_PRIb,beta_p);
        betas(:,i)=beta_it(:);
        fprintf('Replication: %f\n',i)
    end
end    
    
se=std(betas,0,2);

% Alternative p-value calculation. It assumes the distribution under
% the null is the empirical distribution shifted (minus the mean)

% boot_t=(betas-kron(ones(1,nReps),mean(betas,2)))./kron(ones(1,nReps),se);
% coef_t=kron(ones(1,nReps),beta(:)./se);
% pvals=sum(abs(boot_t)>abs(coef_t),2)/nReps;

UCI=prctile(betas,97.5,2);
LCI=prctile(betas,2.5,2);

%p-value computation. It assummes normality of t-stats under null (default in Stata)
pvals=min([2*normcdf(beta_p(:)./se) 2*(1-normcdf(beta_p(:)./se))],[],2);

end